home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / data_smooth / smooth2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  1.7 KB  |  45 lines
  1. /* Program 2. Smoothing and interpolation with second differences. */
  2. /* Contribution to Graphic Gems IV */
  3.  
  4. /* Paul H. C. Eilers, DCMR Milieudienst Rijnmond, 's-Gravelandseweg 565,
  5.    3119 XT Schiedam, The Netherlands, E-Mail: paul@dcmr.nl */
  6.  
  7. #define MMAX 100  /* choose the right length for your application */
  8.  
  9. typedef float vec[MMAX + 1];
  10.  
  11. void smooth2(vec w, vec y, vec z, float lambda, int m)
  12. /* Smoothing and interpolation with second differences.
  13.    Input:  weights (w), data (y): vector from 1 to m.
  14.    Input:  smoothing parameter (lambda), length (m).
  15.    Output: smoothed vector (z): vector from 1 to m. */
  16. {
  17.   int i, i1, i2;
  18.   vec c, d, e;
  19.   d[1] = w[1] + lambda;
  20.   c[1] = -2 * lambda / d[1];
  21.   e[1] = lambda /d[1];
  22.   z[1] = w[1] * y[1];
  23.   d[2] = w[2] + 5 * lambda - d[1] * c[1] *  c[1];
  24.   c[2] = (-4 * lambda - d[1] * c[1] * e[1]) / d[2];
  25.   e[2] = lambda / d[2];
  26.   z[2] = w[2] * y[2] - c[1] * z[1];
  27.   for (i = 3; i < m - 1; i++) {
  28.     i1 = i - 1; i2 = i - 2;
  29.     d[i]= w[i] + 6 * lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
  30.     c[i] = (-4 * lambda -d[i1] * c[i1] * e[i1])/ d[i];
  31.     e[i] = lambda / d[i];
  32.     z[i] = w[i] * y[i] - c[i1] * z[i1] - e[i2] * z[i2];
  33.   };
  34.   i1 = m - 2; i2 = m - 3;
  35.   d[m - 1] = w[m - 1] + 5 * lambda -c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
  36.   c[m - 1] = (-2 * lambda - d[i1] * c[i1] * e[i1]) / d[m - 1];
  37.   z[m - 1] = w[m - 1] * y[m - 1] - c[i1] * z[i1] - e[i2] * z[i2];
  38.   i1 = m - 1; i2 = m - 2;
  39.   d[m] = w[m] + lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
  40.   z[m] = (w[m] * y[m] - c[i1] * z[i1] - e[i2] * z[i2]) / d[m];
  41.   z[m - 1] = z[m - 1] / d[m - 1] - c[m - 1] * z[m];
  42.   for (i = m - 2; 1<= i; i--)
  43.      z[i] = z[i] / d[i] - c[i] * z[i + 1] - e[i] * z[i + 2];
  44. }
  45.